Make three sets of plots for Around the island fDOM presentation
Indices: HIX, A, T, B
Details: Make colors pop (green to red) and drop outliers to get nice spread.
Map
Correlations to ‘distance to shore’
Correlations to silicate, N+N, Phosphorous
library(ggplot2)
library(data.table)
theme_set(theme_minimal())
l_dat = fread("data/clean/ati_fdom_long.csv")
dat = fread("data/clean/ati_fdom_clean.csv")
meta = fread("data/sample_data/islandwide_sites_allMetadata.csv")
inds = c(
'CobleA',
'CobleB',
'CobleT',
'HIX'
)
# log transform all indices
l_dat$log_index_value = log(l_dat$index_value)
# set factor levels for Habitat
Summary of samples collected
site_meta = meta[!is.na(as.numeric(Site)) & May2021 == "Yes"]
## Warning in eval(.massagei(isub), x, ienv): NAs introduced by coercion
site_meta[!(Site %in% dat$Site) & May2021 == "Yes"]
## Site Jan2016 May2016 July2016 July2019 Aug2020 May2021 Latitude Longitude
## 1: 30 Yes Yes Yes <NA> Yes Yes -17.59838 -149.8108
## 2: 97 Yes Yes Yes <NA> Yes Yes -17.54098 -149.9068
## 3: 122 Yes Yes Yes Yes Yes Yes -17.48485 -149.8962
## 4: 125 Yes Yes Yes <NA> Yes Yes -17.49020 -149.8827
## 5: 129 Yes Yes Yes <NA> Yes Yes -17.48973 -149.8727
## 6: 141 Yes Yes <NA> <NA> Yes Yes -17.49316 -149.8674
## 7: 156 Yes Yes <NA> <NA> Yes Yes -17.48503 -149.8313
## 8: 207 <NA> <NA> <NA> <NA> Yes Yes NA NA
## 9: 208 <NA> <NA> <NA> <NA> Yes Yes NA NA
## 10: 15 <NA> <NA> <NA> <NA> <NA> Yes NA NA
## 11: 210 <NA> <NA> <NA> <NA> <NA> Yes NA NA
## Distsance.to.crest Distance_to_shore Distance_to_pass
## 1: 716.16 718.52 3507.52
## 2: 194.05 563.28 160.22
## 3: 162.06 576.72 88.94
## 4: 789.48 137.46 1347.58
## 5: 618.38 232.22 930.83
## 6: 729.33 123.28 589.73
## 7: 843.02 44.33 787.09
## 8: NA NA NA
## 9: NA NA NA
## 10: NA NA NA
## 11: NA NA NA
## Distance_to_deep_lagoon_water Distance_to_population_center
## 1: 389.04 718.52
## 2: 99.38 7334.53
## 3: 274.81 2147.15
## 4: 371.65 592.26
## 5: 65.99 232.22
## 6: 31.14 123.28
## 7: 19.34 2786.91
## 8: NA NA
## 9: NA NA
## 10: NA NA
## 11: NA NA
## Population_center Island_shore Habitat Lagoon
## 1: Ma'atea SW Mid lagoon Maatea
## 2: Vai'anae/Atiha SW Reef pass Haapiti
## 3: Papeto'ai N Reef pass Hauru
## 4: Papeto'ai N Fringing reef Papetoai
## 5: Papeto'ai N Fringing reef Papetoai
## 6: Papeto'ai N Fringing reef Papetoai
## 7: Paopao N Fringing reef Pihaena
## 8: <NA> <NA> <NA> <NA>
## 9: <NA> <NA> <NA> <NA>
## 10: <NA> <NA> <NA> <NA>
## 11: <NA> <NA> <NA> <NA>
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.1; sf_use_s2() is TRUE
library(sp)
library(OpenStreetMap)
library(viridisLite)
map_long_dat = l_dat[!is.na(Latitude)]
# make spacial points data frame
map_sp = SpatialPointsDataFrame(data = map_long_dat,
coords = list(map_long_dat$Latitude,
map_long_dat$Longitude))
# Check geographic range of sampling points
limits = c(
min(map_long_dat$Longitude),
min(map_long_dat$Latitude),
max(map_long_dat$Longitude),
max(map_long_dat$Latitude)
)
# define a bounding box with a small cushion around the minimum and maximum
bbox = list(
xmin = limits[1] - 0.03,
ymin = limits[2] - 0.04,
xmax = limits[3] + 0.03,
ymax = limits[4] + 0.04
)
# get basemap
sa_map <- openmap(c(bbox$ymax, bbox$xmin),
c(bbox$ymin, bbox$xmax),
type = "stamen-terrain",
mergeTiles = TRUE)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
sa_map2 <- openproj(sa_map)
mo_map = function(an_index, outlier_n){
# remove highest and lowest n samples
ind_map_dat = map_long_dat[index_name == an_index &
!(is.infinite(log_index_value))][order(log_index_value)]
ind_map_dat = head(ind_map_dat, n = -(outlier_n))
ind_map_dat = tail(ind_map_dat, n = -(outlier_n))
sa_map2_plt = OpenStreetMap::autoplot.OpenStreetMap(sa_map2)+
geom_point(data = ind_map_dat,
aes(x = Longitude,
y = Latitude,
color = log_index_value))+
labs(title = paste("Log",an_index),
x = "Lon",
y = "Lat")+
scale_color_gradient(low = "green", high = "red")+
guides(fill=guide_legend(title=NULL))
print(sa_map2_plt)
ggsave(filename = paste0("plots/exploration/ATI/map_",an_index,".png"), plot = sa_map2_plt)
}
lapply(inds, mo_map, outlier_n = 10)
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## [[1]]
## [1] "plots/exploration/ATI/map_CobleA.png"
##
## [[2]]
## [1] "plots/exploration/ATI/map_CobleB.png"
##
## [[3]]
## [1] "plots/exploration/ATI/map_CobleT.png"
##
## [[4]]
## [1] "plots/exploration/ATI/map_HIX.png"
# function for making correlation plots
cor_plot = function(ind,
cor_var,
outlier_n,
cor_var_outlier_n){
#subset and drop outliers
ind_dat = l_dat[index_name == ind & !is.infinite(log_index_value)][order(index_value)]
ind_dat = head(ind_dat, n = -(outlier_n))
if(!is.null(cor_var_outlier_n)){
ind_dat = head(ind_dat[order(ind_dat[[cor_var]])], n = -(cor_var_outlier_n))
}
# model
ind_lm = lm(as.formula(paste("log_index_value ~", cor_var)),
data = ind_dat)
ind_sum = summary(ind_lm)
p = ggplot(data = ind_dat, aes_string(x = cor_var,
y = "log_index_value"))+
geom_point(aes(color = Habitat))+
geom_abline(slope = ind_lm$coefficients[2],
intercept = ind_lm$coefficients[1])+
labs(y = paste( "Log",ind),
subtitle = paste("R sq. =",
round(ind_sum$r.squared, 3)))
ggsave(paste0("plots/exploration/ATI/",cor_var,"_",ind,".png"))
return(p)
}
# apply correlation plots for distance to shore
d_shore_ps = lapply(inds,
cor_plot,
cor_var = "Distance_to_shore",
outlier_n = 10,
cor_var_outlier_n = NULL)
## Saving 10 x 5 in image
## Warning: Removed 7 rows containing missing values (geom_point).
## Saving 10 x 5 in image
## Warning: Removed 7 rows containing missing values (geom_point).
## Saving 10 x 5 in image
## Warning: Removed 8 rows containing missing values (geom_point).
## Saving 10 x 5 in image
## Warning: Removed 8 rows containing missing values (geom_point).
d_shore_ps
## [[1]]
## Warning: Removed 7 rows containing missing values (geom_point).
##
## [[2]]
## Warning: Removed 7 rows containing missing values (geom_point).
##
## [[3]]
## Warning: Removed 8 rows containing missing values (geom_point).
##
## [[4]]
## Warning: Removed 8 rows containing missing values (geom_point).
# apply correlation to nutrients
# N+N
NplusN_ps = lapply(inds,
cor_plot,
cor_var = "Nitrite_plus_Nitrate",
outlier_n = 10,
cor_var_outlier_n = 2)
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
NplusN_ps
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
# Phosphate
Phos_ps = lapply(inds,
cor_plot,
cor_var = "Phosphate",
outlier_n = 10,
cor_var_outlier_n = 2)
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
Phos_ps
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
# Silicate
Sil_ps = lapply(inds,
cor_plot,
cor_var = "Silicate",
outlier_n = 10,
cor_var_outlier_n = 4)
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
Sil_ps
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
hab_box = function(ind){
p = ggplot(l_dat[index_name == ind], aes(x = index_name,
y = log_index_value,
color = factor(Habitat, levels = c("Reef pass",
"Reef crest",
"Mid lagoon",
"Bay",
"Fringing reef",
""))))+
geom_boxplot()+
labs(title = paste0("log ",ind))+
facet_wrap(vars(Island_shore), nrow = 2)+
theme(legend.title=element_blank())
print(p)
ggsave(paste0("plots/exploration/ATI/Habitat_boxplot_",ind,".png"))
}
lapply(inds, hab_box)
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Warning: Removed 16 rows containing non-finite values (stat_boxplot).
## Saving 10 x 5 in image
## Warning: Removed 16 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Saving 10 x 5 in image
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## [[1]]
## [1] "plots/exploration/ATI/Habitat_boxplot_CobleA.png"
##
## [[2]]
## [1] "plots/exploration/ATI/Habitat_boxplot_CobleB.png"
##
## [[3]]
## [1] "plots/exploration/ATI/Habitat_boxplot_CobleT.png"
##
## [[4]]
## [1] "plots/exploration/ATI/Habitat_boxplot_HIX.png"